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I .  Introduction 


Due  to  the  current  importance  of  dynamic  analysis  for  civil  engineering 
structures,  especially  with  regard  to  seismic  loading,  there  has  been  a  re¬ 
surgence  of  interest  in  approximate  methods  for  integrating  the  equations  of 
motion.  For  many  years  the  emphasis  has  been  on  the  approximation  of  the 
spatial  deformation  of  these  structures,  using  mor^  sophisticated  finite 
element  models,  and  including  the  effects  of  geometric  and  material  nonlinear¬ 
ity  [I].  Some  of  this  emphasis  is  now  being  aimed  toward  the  understanding  of 
the  characteristics  of  direct  integration  operators. 

In  order  to  be  precise  about  the  scope  of  this  review,  the  distinction 
between  forced  structural  vibration  and  wave  propagation  is  made.  Structural 
vibration  problems  are  almost  always  dominated  by  low-frequency  components  of 
the  response,  since  the  energy  requirements  for  exciting  the  higher  modes  are 
so  severe.  This  is  fortunate  for  the  analyst,  since  the  deterioration  of  the 
accuracy  of  the  frequencies,  and  mode  shapes  makes  the  computed  results  ques¬ 
tionable  for  these  higher  modes.  It  would  seem  reasonable,  therefore,  not  to 
spend  needless  time  and  effort  trying  to  make  the  computed  solution  accurate 
with  respect  to  these  modes. 

High-frequency  response  is  very  important  for  many  wave  propagation 
problems,  however,  especially  where  discontinuities  in  velocity  or  acceleration 
persist.  If  the  response  in  the  region  around  the  discontinuity  is  of  inter¬ 
est,  the  analyst  will  be  forced  to  use  a  time  step  of  integration  so  small  that 
explicit  integration  methods,  such  as  central-differencing  with  a  lumped  mass 
matrix,  become  attractive.  It  should  be  noted  that  wave  propagation  problems 
without  velocity  or  acceleration  discontinuities  are  handled  quite  nicely  by 
the  same  techniques  that  are  used  for  structural  vibration  [2], 
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In  this  review,  several  alternative  methods  for  carrying  out  the  step- 
by-step  integration  of  the  equations  of  motion  of  a  structural  system  will  be 
discussed.  In  the  next  section,  the  characteristics  of  general  structural 
systems  are  emphasized  so  that  the  criteria  for  choosing  a  particular  integra¬ 
tion  operator  can  be  understood  on  these  terms >  Then,  three  of  the  most  popu¬ 
lar  methods  are  analyzed  and  compared.  Following  that  is  a  discussion  of  the 
relationship  between  these  popular  methods  and  a  class  of  methods  referred  to  in 
the  literature  as  stiffly-stable  [3].  Finally,  recent  results  on  the  applica¬ 
tion  of  these  results  to  nonlinear  problems  are  summarized. 

II.  Structural  System  Characteristics 

The  governing  equation  of  interest  is  given  by 

M  (l(t)  +  C  u(t)  +  K  u(t )  =  F(t)  ,  (II. 1) 

*«•  «•  *■»  *»  *» 

•m  ^ 

where  M  ,  C  ,  and  K  are  the  mass,  damping,  and  stiffness  matrices,  respec¬ 
tively;  U(t)  ,  u(t)  ,  and  u(t)  are  the  acceleration,  velocity,  and  displace¬ 
ment  vectors  at  time  t  ,  respectively;  and  F(t)  is  the  force  vector.  An 
incremental  form  of  this  equation,  valid  for  elastic-plastic  constitutive  be¬ 
havior  and  geometrically  nonlinear  behavior,  has  been  developed  in  [4],  This 
incremental  form  is  used,  in  conjunction  with  "residual  load  correction"  in 
order  to  solve  beam  and  axisymmetric  shell  dynamic  problems.  The  governing 
equation  (II.l)  can  be  derived  in  a  number  of  ways,  such  as  by  Rayleigh-Ritz- 
Galerkin  or  other  finite  element  methods,  but  the  mass  matrix  will  always  be 
assumed  to  be  positive  definite. 

For  most  structural  problems,  the  mass,  damping,  and  stiffness  matrices 
* 

are  sparse,  banded,  and  symmetric.  Attempts  to  solve  ill.l),  whether  by  direct 
integration  (step-by-step)  or  by  modal  superposition,  should  take  advantage  of 
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these  characteristics.  The  initial  conditions  for  (II.l)  involve  the  displace¬ 
ment  and  velocity  at  t  =  0  ;  therefore,  the  direct  integration  procedure  re¬ 
quires  at  least  this  amount  of  information  about  the  previous  time,  t  ,  in 
order  to  predict  the  state  of  motion  at  the  current  time,  tR+^  .  Additional 
information  about  the  state  at  time  t  ,  such  as  acceleration,  o’*  about  the 
state  of  motion  at  an  even  earlier  time  (say  t  leads  to:  (a)  special  start¬ 
ing  procedures  for  the  integration  method;  (b)  extra  storage  requirements;  and 
(c)  extraneous  solutions  [5]  that  may  create  accuracy  or  instability  problems. 
Unless  ample  justification  for  the  extra  information  is  given,  the  minimum  amount 
should  be  used.  ‘To  the  writer’s  knowledge,  only  one  such  method  has  been  pro¬ 
posed  [6],  the  Gurtin  Averaging  operator,  that  is  unconditionally  stable;  how¬ 
ever,  the  accuracy  of  the  method  is  low  unless  small  time  steps  or  negative 
damping  is  used. 

The  characteristics  of  (II.l)  are  best  defined  in  terms  of  its  natural 
frequencies  and  mode  shapes.  For  a  continuous  system,  the  frequency  spectrum 
ranges  from  the  lowest,  or  fundamental,  frequency  up  to  an  infinite  limit  point 
(there  may  be  zero  frequencies  if  rigid-body  modes  are  present).  For  a  dis¬ 
cretized  system,  the  infinite  limit  point  no  longer  exists;  instead,  a  frequency 
exists  that  corresponds  to  the  most  rapidly  varying  (in  space)  mode  shape.  This 
frequency  is  called  the  cut-off  frequency.  If  the  discretized  system  is  excited 
by  forcing  functions  having  frequency  content  above  the  cut-off  frequency,  such 
as  might  be  induced  in  a  wave  propagation  problem,  noise  (random  spatial  response) 
is  generated  in  the  cut-off  modal  response. 

Corresponding  to  each  possible  mode  of  vibration  and  frequency  of  the 
structural  system  is  a  period  of  vibration,  or  time  constant,  and  the  range 
of  periodicity  is  a  measure  of  the  "stiffness”  [3]  of  the  system.  The  largest 


period  is  associated  with  the  fundamental  mode  and  the  smallest  is  associated 
with  the  cut-off  frequency  or,  therefore,  the  "stiffest"  system  component. 

As  Jensen  has  pointed  out  [7],  the  integration  of  (II.l),  which  is,  in 
general,  a  stiff  system,  presents  the  analyst  with  a  dilemma.  If  the  time  step 
is  reduced  in  order  to  accurately  integrate  the  stiff  components,  then  the  step 
will  be  much  too  small  for  the  longer  period  (lower  frequency)  responses, 
resulting  in  excessive  computer  time  for  the  calculations.  On  the  other  hand, 
if  the  time  step  is  chosen  with  regard  for  the  low-frequency  response,  insta¬ 
bility  will  result  for  explicit  methods  and  integration  error  will  give  inac¬ 
curate  solutions  for  the  stiff  components. 

Four  solutions  of  interest  can  be  identified  here:  (a)  the  exact  solu¬ 
tion  of  the  equations  of  motion  of  the  continuous  system;  (b)  the  exact  solu¬ 
tion  of  (II.l),  using  infinite-precision  arithmetic;  (c)  the  exact  solution  of 

(11.1) ,  using  finite-precision  arithmetic;  and  (d)  the  approximate  solution  of 

(11.1) ,  using  finite-precision  arithmetic  and  a  direct  integration  operator. 
Truncation  error,  such  as  the  noise  referred  to  above,  represents  error  incurred 
from  the  truncation  of  a  continuous  to  a  discretized  system  and  is  the  differ¬ 
ence  between  (a)  and  (b).  Round-off  error  depends  on  the  degree  of  precision 
of  the  computational  arithmetic  and  is  the  difference  between  (b)  and  (c). 
Integration  error  is  the  difference  between  (c)  and  (d). 

The  primary  motive  for  studying  various  direct  integration  operators 
is  to  understand  the  nature  of  integration  error,  how  to  control  it,  and  how 
to  estimate  the  computing  costs  of  this  control.  An  integration  operator  [8] 
is  defined  as  a  transformation  on  the  acceleration,  velocity,  and  displacement 
vectors  at  time  t^  (and,  possibly,  at  earlier  times)  to  the  acceleration, 
velocity,  and  the  displacement  vectors  at  time  t^+^  .  If  the  mass,  damping, 


and  stiffness  matrices  do  not  depend  on  the  displacements,  or  their  space  and 
time  derivatives,  the  transformation  is  said  to  be  linear.  If  the  transforma¬ 
tion  does  not  depend  on  the  state  at  times  earlier  than  t  ,  it  is  called  a 
single-step  method;  otherwise,  it  is  a  multi-step  method. 

Direct  integration  operators  are  generally  written  in  the  form 

K.  u(t  .)  =  F  +  K  u(t,t  -,...)  ,  (II. 2) 

„  n+i  -  _o  ,  n  n-i 

mt 

where  K,  and  K  are  the  matrices  that  define  the  transformation,  F  is  a 

.1  _o  . 

vector  of  forcing  functions,  and 
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(II. 3a) 
(II. 3b) 


The  superscript  T  indicates  the  transpose  of  a  vector  or  a  matrix.  The 
matrices  and  Kq  are,  in  general,  functions  of  the  mass,  damping  and  stiff¬ 
ness  matrices,  as  well  as  the  time  step  size,  At  =  t^+1  -  t^  .  If  the  matrix 
K.  can  be  put  into  upper  or  lower  triangular  form,  the  operator  is  said  to  be 
explicit;  otherwise,  it  is  implicit.  The  amplification  matrix  of  the  operator 
is  defined  by 

A  =  K ._1  K  ,  (II. 4) 

.  .1  .o 

assuming  that  an  inverse  exists.  Then 

u(tn+1)  =  G  +  Au(tn,tn_  . . )  ,  (II. 5) 


where 


ways.  Lax  and  Richtmyer  [8]  discuss  the  spectral  radius  of  the  amplification 
matrix;  for  structural  dynamics,  this  spectral  radius  is  defined  in  terms  of 


the  resonant  structural  frequencies  of  the  system  and  the  integration  operator 
is  stable  for  all  time  step  sizes  that  cause  the  spectral  radius,  to  be  bounded 
by  unity.  Dahlquist  [9]  has  taken  a  different  approach,  introducing  the  con¬ 
cept  of  A-stability  for  systems  of  first-order  ordinary  differential  equations. 
By  this  he  means  that  the  error  introduced  into  the  approximate  solution  by  the 
particular  integration  method  remains  uniformly  bounded  for  any  time  step  size. 
This  coincides  with  the  Lax-Richtmyer  notion  of  unconditional  stability;  i.e., 
the  spectral  rtdius  of  the  operator  is  bounded  by  unity  for  all  choices  of 
time  step  size.  Dahlquist  has  proven  that  linear  multi-step  methods  are 
implicit  if  they  are  A-stable  and  that  the  "trapezoidal"  rule  has  the  smallest 
asymptotic  error  of  any  A-stable  method. 


III.  Operator  Comparisons 

In  this  section  the  characteristics  of  the  most  popular  direct  integra¬ 
tion  operators  are  discussed  and  compared.  The  central  difference  operator, 
for  example,  has  been  shown  to  be  conditionally  stable  [10]  and,  in  addition, 
Krieg  [11]  has  found  that  no  explicit  integration  operator  of  order  two  has  a 
stability  region  greater  than  the  central,  difference  operator.  The  Houbolt 
operator  [12,13],  on  tne  other  hand,  is  unconditionally  stable  [14]  and  has 
been  compared  to  the  central  difference  curator  for  accuracy  and  speed  [15]. 
Both  the  Newmark  [16]  and  Wilson  Averaging  [17,18]  operators  have  alt  been 
shown  to  be  unconditionally  stable  [6]  and  comparisons  ha\e  been  made  wi-*-h  a 
precise  integration  operator  [19]  based  on  modal  superposition.  In  spit?  of 
this  work,  however,  some  direct  comparisons  seem  to  be  in  order. 


'i 

3 

y 


i 

* 

i 


1 

I 

I 

it 

I 


| 


i 


a 


-  7  - 


The  Houbolt  operator  is  obtained  by  fitting  a  cubic  polynomial  through 
values  of  the  current  displacement  (to  be  found)  and  the  three  previous  values. 
This  necessitates  a  special  starting  procedure.  Substituting  into  (II.l)  the 
first  and  second  derivatives  of  this  polynomial,  evaluated  at  the  current  time, 
results  in  an  equation  of  motion 

<?  +  §  +  7  “^1  * 

+  (|  M  +  4  AtC)u  -  (2M  +  |  AtC)u  ,  (III.l) 

2.2.. n  .4.  .n-l 

•«*  • 

.  (i  «  .  i  4tC)Vj  . 

•»  <v 

Following  the  s  U'oility  procedure  of  von  Neumann  [20],  let 

un  =  Xn  d  ,  (III. 2) 

where  d  is  an  arbitrary  error  vector.  Then,  the  characteristic  equation 
for  A  is 

U  +  |  C  +  ~  n)A3  -  (f  +  |  n)A2  +  (2  +  |  n)X  -  (|  +  i  n)  «  0  ,  (III. 3) 

2  2  '1 

where  g  s  u  (At)  ,  n=  AtM  C  ,  to  is  any  of  the  undamped  structural 
resonant  frequencies,  and  At  is  the  time  step  size.  The  characteristics  of 
the  three  roots  of  this  equation  (one  real,  two  complex  conjugates)  can  be 
investigated  in  terms  of  a  single  degree- of- freedom  resonance  and  a  damping 
factor  in  order  to  determine  the  extent  of  artificial  damping  and  periodicity 
error.  Figure  1  shows  the  values  of  the  moduli  of  these  roots,  for  the  case  of 
zero  damping,  plotted  against  the  variable  to  (At)  .  The  extraneous  root  [5], 
R  ,  causes  the  greater  artificial  damping  (no  artificial  damping  corresponds 
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to  =  1)  but  it  seems  clear  that  at  least  fifteen  time  steps  (Levy  and 
Kroll  [13]  suggest  thirty)  are  required  to  maintain  the  modal  amplitude  at  any¬ 
where  near  the  correct  value.  Figure  2  shows  the  same  roots,  plotted  against 

2  2 

smaller  values  of  u  (At)  ,  showing  that  the  extraneous  root  becomes  the  major 
factor  in  the  artificial  damping  of  the  lower  modes  of  structural  response. 

As  an  example  of  the  choice  of  time  step  size  to  be  used  in  conjunc¬ 
tion  with  the  Houbolt  operator,  consider  the  spherical  shell  cap  under  a  point 
load  at  the  apex  as  analyzed  by  McNamara  and  Marcal  [4],  The  load  is  applied 
suddenly  and  maintained  at  a  constant  value  throughout  the  analysis.  Although 
the  shell  is  so  thin  and  shallow  that  geometric  nonlinearity  dominates  the 
dynamic  response,  some  insight  into  the  effect  of  artificial  damping  and  peri¬ 
odicity  error  can  be  gained  by  comparing  the  time  step  used  in  the  nonlinear 
-5 

analysis  (At  =  10  seconds)  with  the  periods  of  vibration  for  the  lowest 
(linear)  modes.  Such  a  comparison  is  given  in  Table  I.  Noting  that  the 
abscissa  of  Fig.  2  can  be  written  as  (At/2mT)  ,  where  T  is  the  period  of 


TABLE  I 

Spherical  Shell  Cap  Modes  and  Periods 


Mode 

Frequency,  Hz 

Period,  usee 

At/Period 

1 

892.2 

1122.0 

.00892 

2 

1165.0 

858.  C 

.01165 

3 

.1853.0 

540.0 

.01853 

4 

3092.0 

323.5 

.03092 

,  the 

error  due  to  artificial  damping  after  1000 

time  steps  : 

fourth  mode  response  is  about  two  percent.  This  can  be  considered  negligible. 
It  is  unlikely  that  modal  response  above  the  fourth  mode  will  be  adequately 
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represented,  however,  unless  the  geometric  nonlinearities  increase  the  periods 
significantly < 

The  Wilson  Averaging  operator  has  also  been  analyzed  [6]  for  artificial 
damping  and  comparisons  with  the  Houbolt  operator  are  shown  in  Figs.  3  and  4. 
While  ther'  is  still  strong  damping  for  higher  modes,  the  increased  accuracy 
of  the  Wilson  Averaging  operator  is  evident.  It  should  be  noted,  at  this 
juncture,  that  the  Newmark  generalized  acceleration  operator,  with  y  =  1/2  , 

(a)  has  a  zero  extraneous  root  and  (b)  has  no  artificial  damping,  regardless 
of  the  value  of  3  [6]. 

The  periodicity  error  comparisons  are  given  in  Figs.  5  and  6.  The 
bust  performance  is  again  exhibited  by  the  Newmark  operator  with  3  =  1/4 
(note  that  3=0  causes  the  computed  periods  to  be  smaller  than  the  exact 
periods,  in  contrast  to  the  other  integration  operators).  Of  the  two  operators 
with  nonvanishing  extraneous  roots,  the  Wilson  Averaging  operator  is  again 
superior.  Figure  7  depicts  the  influence  that  real  damping  has  on  the  artifi¬ 
cial  damping  in  the  Houbolt  operator.  The  effect  on  the  complex  conjugate  roots 
is  mild,  since  linearity  of  the  change  in  damping  is  indicated  for  reasonable 
values  of  true  damping.  More  disturbing  is  the  effect  displayed  by  the  extrane¬ 
ous  root,  which  shows  a  nonlinear  decrease  in  damping  as  a  function  of  increased 
true  damping,  even  indicating  an  instability  for  large  enough  true  damping. 
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IV.  Stiffly  Stable  Methods 

Gear  [3]  has  suggested  that  the  requirement  for  absolute  stability 
for  all  components  of  the  solution,  regardless  of  the  time  constant  of  the 
component,  is  too  restrictive  and  has  indicated  a  preference  for  "stiffly 
stable"  methods,  i.e.,  methods  which  have  regions  of  stability  sufficient  to 
include  frequencies  up  to  the  cut-off  frequency.  From  Fig.  8,  the  character¬ 
istics  of  a  stiffly  stable  method  are:  (a)  stable  and  accurate  solutions  in 
the  cross-hatched  region  A;  (b)  stable,  but  not  necessarily  accurate,  solutions 
in  the  cross-hatched  region  B;  and  (c)  solutions  of  questionable  stability  and 
accuracy  elsewhere.  Each  particular  stiffly  stable  method  has  it'i  own  generic 
parameters  £  ,  n  and  £  that  define  the  regions  of  relative  accuracy  and 
stability.  The  A-stable  methods  correspond  to  n  =  0  . 

The  argument  for  using  stiffly  stable  methods  has  been  made,  for  struc¬ 
tural  dynamics  problems,  by  Jensen  [73.  He  suggests  that  the  order  of  the 
stiffly  stable  method  be.  varied,  depending  upon  the  particular  problem  being 
solved,  in  order  to  maintain  accuracy  for  the  important  components  of  the 
solution.  It  should  be  noted  that  stiffly  stable  methods  of  high  order  have 
the  disadvantage  of  implicit  backward  difference  operators  such  as  the  Houbolt 
operator — namely,  that  considerable  storage  space  is  occupied  by  the  vectors 
of  past  displacement ,  velocity  or  acceleration. 

Linear  multistep  methods  for  a  system  of  first-order  ordinary  differen¬ 
tial  equations 

•  _ 
u  = 

with  initial  conditions 

u(tQ) 

f 

1 


f(u,t) 


=  u 

-O 


(IV. 1) 


(IV. 2) 
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are  written  in  the  form  [21] 

Jj  “i  ?n.i  =  4t  j/i  4n+1  ,  n  =  -1  ,  0  ,  1  ....  (IV. 3) 

where  At  is  the  step  size,  u  =  u(t  +  nAt)  ,  and  o  i  0  .  If  B  =  0  , 

-n  -  o  m  m 

(IV. 3)  is  an  explicit  operator  for  u  .  and  is  referred  to  as  a  predictor;  if 

_n+m 

/  0  ,  (IV. 3)  is  an  implicit  operator  and  is  referred  to  as  a  corrector. 
Stiffly  stable  methods  fall  into  the  latter  class;  therefore,  it  is  of  interest 
to  find  out  how  they  compare  with  implicit,  unconditionally  stable  operators, 
such  as  the  Newmark,  Houbolt  and  Wilson  Averaging  operators. 

In  order  to  make  this  comparison,  the  equations  of  motion,  (II.l),  must 
be  transformed  from  n  second-order  ordinary  differential  equations  into  2n 
first-order  ordinary  differential  equations.  Jensen  [7]  has  given  such  a  deconv 
position  that  has  the  advantage  of  retaining  the  same  size  for  the  system  of 
implicit  equations,  i.e.,  the  additional  n  equations  are  already  in  upper  or 
lower  trianguiarized  form. 

Thus,  define 

v(t)  =  M  u(t)  +  (C  +  yAtK)u(t)  ,  (IV. 4) 

where  Y  =  B  /a  .  Then 
m  m 

v(t)  =  F(t)  -  Ku(t)  +  yAtKu(t)  .  (iv.s) 


As  an  example,  suppose  that  a  single-step  (m=2)  method  is  being  used.  Then 
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;]  fO  (O  (O 

"v"  .  =  °h-  +  BAtrr  +  YAtU-J 

C  -Jn+1  -  Jn  C  -  Jn  (  -Jn 


(IV. 7) 


The  governing  system  to  be  solved  consists  of  equations  (IV.4),  (IV. 5) 
and  (IV.7).  Some  care  must  be  taken  when  evaluating  the  stability  properties 
"this  operator,  however,  m  that  the  starting  procedure  must  be  taken  into 
account.  Since  u(tQ)  and  u(tQ)  are  the  prescribed  initial  conditions,  the 
initial  values  for  v(to)  and  v(tQ)  can  be  found  from  (IV.4)  and  (IV.5): 


v(t  )  =  Mu(t  )  +  (C  +  yitK)u(t  ) 


(IV. 8) 


v(to)  =  F(to}  “  Ku(to>  +  YAtKu(t  )  .  (IV. 9) 

%  v  •*  v  ***  v  %%  O 


Therefore,  the  second  of  the  partitioned  equations  (IV.7)  must  be  rewritten, 
in  general. 


+  (oH  +  6r*t2K)u 
-  -  -n 


(IV. 10) 


+  (aC  +  (ay-S)AtK)u 
-  «,  «.n 


The  governing  system  to  be  solved  can  then  be  written  in  the  Lax- 
Richtmyer  form  (see  (II. 2))  for  direct  integration  operators  where 
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and 


[^3  « 


0  -yAtl 


oy  ■ 


ctl 


-yAtl 


-C-yAtK  I  -M 


0  -yAtK 


BAtl 


aC+(ay-e)AtK  0  aM+ByAt  K  0 


(DT  .  <0,  B4tFn ,  0,  Fn+1  , 


f  =  ^u,  v,  o,  $  y 


(IV.ll) 


(IV. 12) 


(IV. 13) 


where  I  is  the  identity  matrix.  Note  that 

\ 

(IV. 14) 

Without  too  much  difficulty  the  inverse  of  (IV.ll)  can  be  found  and 


premultiplied  by  (IV.12).  Then  the  amplification  matrix  is  given  by 
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D”1(aH+OYAtC-$YAt2K)  0  D“1(aY+B)AtM 


[A]  = 


aC-0AtK 


D^crr+BMtK 


0  D“1(aH-eAtC-6Y<it2K)  0 


-BAtK 


,  (IV. 15) 


where  D  =  M+yAtC+y  At  K  .  This  expression  confirms  the  fact  that,  if  y  =  0 
and  the  mass  matrix  is  diagonal,  the  procedure  is  explicit.  With  y  =  o  and 
a  distributed  mass  matrix,  the  operator  remains  implicit,  since  the  mass  matrix 
must  be  inverted. 

The  spectral  character  of  this  operator  can  be  investigated  by  finding 
the  eigenvalues  of  the  amplification  matrix  as  function  of  the  structural 
frequencies. 

The  characteristic  equation  for  A  is 

\  1/ II  •  i  1  M  A.  I  1  2,.\  I  *1  I  . .  A  .  .  .  . 


X  [D  (aM+ayAtC-ByAt  K)  -  X]  [D  (aM-BAtC-BYAt  K)  -  X] 


+  X2D"1KD“1H(aY+0)2At2  =  0  . 


(IV. 16) 


IWo  of  the  eigenvalues  are  therefore  seen  to  vanish,  indicating  the  presence  of 
extraneous  roots  [5]  in  the  operator,  as  might  have  been  expected.  For  the  case 
of  zero  damping,  the  remaining  two  roots  are  complex  conjugates,  given  by 


Xlj2  =  D_1(aM-8YAt2K  +  iAt(aY+0)/iTM  )  , 


(IV. 17) 


or,  in  terms  of  the  structural  frequencies. 
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(IV. 18) 


It  would  be  interesting  to  compare  this  result,  specialized  for  the  trapezoidal 
rule  (b  :  1,  y  *  M  1/2)  ,  with  the  conventional  integration  operators  of 
structural  dynamics.  In  this  case, 


,  1  2.  .2  ,  .  . . 

1  -  jj'  u  At  ±  iwAt 

1  5  2 

1  +  i  w  At 


(IV. 19) 


This  result  turns  out  to  be  identical  [6]  with  the  result  for  the  Newmark  gener¬ 
alized  acceleration  method  with  y  =  1/2,  8  =  1/4  .  Since  the  trapezoidal  rule 
was  shown  by  Dahlquist  [9]  to  have  the  smallest  asymptotic  error  of  all  order 
two  methods  of  this  class t  in  addition  to  being  A-stable,  it  seems  unlikely 
that  much  improvement  can  be  made  with  other  values  of  the  parameters. 

V.  Nonlinear  Problems 

Other  stiffly  stable  operators  of  higher  order  can  be  formulated 
[22-24],  but  since  the  storage  requirements  are  large  and  since  these  higher 
order  methods  are  also  implicit,  there  seems  to  be  little  motivation  for  their 
study,  unless  it  can  be  shown  that  the  error  in  the  Newmark  method  is  excessive. 
Stricklin,  et  al.  [25]  have  indicated  a  more  serious  problem — that  the  Newmark 
method  degenerates  when  nonlinear  problems  are  being  analyzed,  leading  to  un¬ 
stable  solutions.  Both  [25]  and  [4]  have  adopted  the  Houbolt  method  for  this 
class  of  problems  in  order  to  take  advantage  of  the  artificial  damping  of 
spurious  components  in  the  solution  due  to  the  transpostion  of  nonlinear  terms 
to  the  right-hand  side  of  (II. 1).  Others  [26,27,19]  have  suggested  the  use  of 
either  limited  or  unlimited  modal  superposition  methods,  even  for  nonlinear 


itnaxow 
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problems.  (It  is  worth  noting  that  the  Newmark  operator  has  been  used  very 
successfully  for  stress  formulations  of  the  equations  of  motion  [28].) 

Boggs  [29]  has  recently  shown  that  the  trapezoidal  rule  (or  Newmark' s 
method)  is  an  effective  procedure  for  solving  nonlinear  equations,  provided  that 
proper  predictor-evaluation-corrector  algorithms  (PEC)  are  chosen.  This  means 
that  a  solution  is  predicted  on  the  basis  of  stiffness  matrices  (initial  stress, 
initial  displacement,  and  small  displacement  matrices  [30])  evaluated  at  the 
end  of  the  last  step;  then,  these  matrices  are  re-evaluated  on  the  basis  of  the 
predicted  solution;  finally,  a  corrected  solution  is  sought.  Boggs  has  found 
that  an  explicit  predictor,  an  evaluation  based  on  the  predicted  solution,  and 
the  trapezoidal  rule  corrector  proves  to  be  adequate.  He  also  explored  iterative 
methods  which  avoided  the  inversion  of  the  Jacobian  (iterative  explicit). 

In  another  recent  publication.  Weeks  [31]  has  evaluated  both  the  trape¬ 
zoidal  rule  and  the  Houbolt  operator  (as  well  as  central  differencing)  for 
geometrically  nonlinear  dynamic  structural  response  problems.  He  found  that  a 
Newton  Raphson  iterative  technique  for  both  operators  led  to  adequate  results, 
and  that  the  pseudo- load  extrapolation  procedure  [25]  caused  instabilities  for 
unconditionally  stable  operators  when  larger  time  steps  are  used.  The  extra 
storage  and  cost  associated  with  the  Newton  Raphson  technique  makes  its  use  in 
practical  situations  doubtful. 

From  all  available  evidence,  then,  it  would  seem  that  the  trapezoidal, 
or  Newmark,  operator  is  the  most  attractive  direct  integration  operator  for 
both  linear  and  nonlinear  problems.  The  suggestions  of  Boggs  [29]  have  very 
nearly  been  applied  by  McNamara  and  Marcal  [4],  who  use  an  implicit  predictor 
(in  this  case,  the  Houbolt  method),  evaluate  on  the  basis  of  the  predicted 
solution,  then  apply  a  "load  correction"  for  the  next  implicit  prediction, 
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based  on  the  residual  error  from  the  equations  of  motion.  This  same  procedure 
should  also  be  applicable  in  conjunction  with  the  trapezoidal  rule  and  would 
seem  to  represent  the  optimum  choice  of  a  direct  integration  operator,  consid¬ 
ering  economy ,  accuracy  and  stability. 
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FIG.  5  PERIODICITY  ERROR  COMPARISON  FOR 
NEWMARK  AND  HOUBOLT  OPERATORS 
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FIG.  7  EFFECT  OF  TRUE  DAMPING  ON  ARTIFICIAL 
DAMPING,  HOUBOLT  OPERATOR. 


